Non-equilibrium quantum domain reconfiguration dynamics in a two-dimensional electronic crystal and a quantum annealer

Relaxation dynamics of complex many-body quantum systems trapped into metastable states is a very active field of research from both the theoretical and experimental point of view with implications in a wide array of topics from macroscopic quantum tunnelling and nucleosynthesis to non-equilibrium superconductivity and energy-efficient memory devices. In this work, we investigate quantum domain reconfiguration dynamics in the electronic superlattice of a quantum material using time-resolved scanning tunneling microscopy and unveil a crossover from temperature to noisy quantum fluctuation dominated dynamics. The process is modeled using a programmable superconducting quantum annealer in which qubit interconnections correspond directly to the microscopic interactions between electrons in the quantum material. Crucially, the dynamics of both the experiment and quantum simulation is driven by spectrally similar pink noise. We find that the simulations reproduce the emergent time evolution and temperature dependence of the experimentally observed electronic domain dynamics.

1 Supplementary Notes 1 1.1 The experimental setup, domain creation by laser or charge injection through tip Domains may be created at 4 K either by exposure to a laser pulse 1 , or by an electrical pulse from an STM tip [2][3][4] , or by a current pulse passed through a set of electrodes on the surface 4,5 .In this work the domain state was created by using an electrical pulse from the STM tip on a freshly cleaved 1T-TaS2 sample in UHV.The tip is initially retracted at 50 nm from the sample's surface, then approached with a bias voltage of 5 V and a tunneling current of 10-50 nA.Alternatively, the tip can be approached to the surface with a voltage bias of 0.05 V and a tunneling current of 1-3 nA, then a 4 V electrical pulse of 100 ms duration is applied.The two described procedures generate a domain state area of around 10 4 nm 2 .The properties of the resulting domain formation 2,3,[6][7][8] , ordering 9 and resistance relaxation as a function of temperature and on different substrates have been reported in detail elsewhere 5 .

Calculation of the STM tip and sample temperature during transformation and scanning
Thermal simulations were performed by solving partial differential equations governing heat conduction using the finite element method (FEM).The thermal model was coupled to the electromagnetic model which calculated the resistive losses (Joule heating) in the contact region:     − ∇ • (∇) =   where   = |∇| 2 .The density , heat capacity   10 thermal conductivity  11 , and electrical conductivity  of all the materials were taken from literature and/or calibrated by experiment 1 .The exact geometry of the W tip, tip mount, tip mount heat sink, 1T-TaS2 crystal sample and sample holder mount in the STM UHV chamber were reproduced in the model calculations.The 1T-TaS2 sample (100  thick) is glued to the conducting (Mo) sample holder substrate, which is clamped to the cryostat vessel at base temperature (4.5 K) with conducting silver paste.The calculation of temperature of the W tip and 1T-TaS2 sample for different tunneling currents and tip penetration depths is performed using the COMSOL package.The geometric elements at the tip-sample contact used in the calculation are shown in Supplementary Fig. 1.
The calculations take into account the temperature of the sample and the tip when a current is passing between the tip and the substrate.The model calculation also takes into account: the temperature dependence of thermal conductivity, and the specific heat capacity of both 1T-TaS2 and tip material (W) the temperature-dependence and anisotropy of electrical conductivity of both materials.For 1T-TaS2, the different electrical conductivities of the C and the H phase are taken into account.thermal clamping of the tungsten rod at a distance 2 mm to the STM tip mount the heat sink of the piezoelectric STM actuator which is in wide-area thermal contact with the tip mount is at base temperature, connected to the cryostat cold finger by multiple copper braids.The thermal contact between the sample and the sample holder is assumed to be ideal.The calculations are shown in Supplementary Fig. 3 for various penetration depths p.As expected, the behavior beyond a few nm is not very dependent on p.A calculation of the tip and sample temperatures in the contact area during single pulse excitation is shown in Supplementary Fig. 2. The maximum temperature at the center reaches ~ 30 K. The temperature drops rapidly away from the contact point.The thermal perturbation is therefore too small to heat the sample through the nearest phase transition to the trigonal (T) or nearly-commensurate (NC) phase.We can safely conclude that the main effect is in response to charge injection, and not thermal.
Supplementary Figure 3 | Vertical (a), and lateral (b), temperature profiles within the 1T-TaS 2 sample for different STM tip penetration depths p (see Supplementary Figs.1,2) resulting from the pulse causing the transformation to the domain state (5V, 50nA).A steady state is assumed, so these are maximum attainable temperatures during the passing of the applied current.
1.2.2Heating of the sample by the tip current during relaxation measurements.
The temperature calculated for the current and voltage parameters used in the relaxation measurements shown in Fig. 1 of the main text are shown in Supplementary Fig. 4. The maximum calculated temperature rise is < 0.000001 .We conclude that the temperature rise due to the STM tip currents used in relaxation measurements can be neglected.
Supplementary Figure 4 | Vertical (a), and lateral (b), temperature profiles within the 1T-TaS2 sample for different STM tip penetration depths p (see Supplementary Figs. 1, 2) during relaxation measurements using the parameters presented in Fig. 1 of the main text.The maximum temperature rise is < 0.000001 K.The sample and tip mount base temperatures are both 4.5 K. .Integrating   (, ) over energy in the interval ±Δ/2, where Δ is the full bandwidth of state   , we obtain a real-space electron density () for state   .The total charge  in state   is given by the integral over a contour  around a localized state ∫ (r)r =   .If charge is localised, maxima in STM images can be interpreted as an electron count, which can be used to quantify the domain wall reconfiguration.By measuring (r) patterns as a function of time, one obtains a quantitative measure of charge reconfiguration in each frame.We note that this does not imply any single particle motion between frames, but measures the number of sites involved in the domain reconfiguration in a given time interval.Furthermore, the absolute number density is not of interest, only the changes with time.The criterion for domain wall characterization is thus somewhat arbitrary.To obtain an accurate rate it is important that the counting is performed consistently.
In these measurements, the border between the ground state and the domain state is first found.Supplementary Fig. 5 shows an example of such an STM image.An area is typically chosen for analysis with as few imperfections and structural defects (such as the bright yellow spots in Supplementary Fig. 5) as possible.A series of N images of the same area (40x40 nm 2 ) are taken at regular intervals to record the relaxation of the domain structure (Supplementary Fig. 6).Tipscanning parameters are -0.8V for voltage bias, 50 pA for tunneling current and 84 nm/s is the tip speed.The tip scans each line from left to right along the x-axis and then moves to the upper line (y-axis).8 minutes are required to scan the entire area under these conditions.Images within the same series are rigidly translated referring to a fixed point common to each image (e.g. a stationary polaron) and then uneven borders are cut.This corrects possible shifts of the tipsample relative position.To compare two consecutive images (i and i+1; with i = 1, . .., N), we detect the polarons' centers in each image (shown by red dots in Supplementary Fig. 7) by using a blob detection algorithm in Python environment.Then the centers' patterns are juxtaposed (Supplementary Fig. 8a) to count the number of polarons which have moved.
The allowed displacement of the electrons is defined by the internal structure of 1T-TaS2 polarons, each of which involves thirteen Ta atoms arranged in the shape of a David's star, with 6 Ta atoms on the vertexes, 6 at the corners, and one in the center (Supplementary Fig. 8b).Due to the periodicity of the system, there are only two non-trivial ways to move a polaron, i.e. move the central atom to a corner or to a vertex of the David's star.This corresponds to a minimum polaron's movement of ~0.33 nm.Thus, to determine if a certain polaron moved between two images, we take its coordinates in the two images and we evaluate if they differ by more than 0.33 nm.In the affirmative case, the polaron is considered as moved.Supplementary Fig. 8c shows one possible arrangement of the possible polaron positions at the domain wall enumerated by Ma et al 2 .The movement of a domain is associated with an energy cost Δ.To first order, ~, where  is the number of polarons involved in the reconfiguration.The fraction (f) of moved polarons is calculated as the number of moved polarons divided by the total number of polarons: . ( We estimated the error of f as given mainly by polarons which are detected only in one of the two compared images (this usually occurs on the border of the images or inside the domain walls).Its value is calculated as the relative difference between the total detected polarons in the two images: The described procedure is repeated for each pairs of consecutive images in a certain series.
Finally, summing all the calculated fractions and dividing by the time (T) needed to scan the N images of the series ( =  × 8 .),we obtain the rate R of movements in the series:

𝛿𝑅 = Σ 𝑖=1…𝑁−1 𝛿𝑓(𝑖, 𝑖 + 1) 𝑇
With this definition, R is always ≤ 1.This is due to the fact that if a polaron moves two (or more) times between two images, it would still be detected as a single event (or 0, if the polaron returns to its starting site).Note that the absolute value of N is of no significance, while the accuracy of the rate measurement is ensured by the use of a consistent algorithm throughout the analysis.
Supplementary Fig. 9 shows the raw data on relaxation at different temperatures.The scatter of the data is a result of the fact that domain reconfiguration process is measured on a relatively small area which is comparable to the domain size.Measurement of a larger area would require faster scanning, which in turn would sacrifice resolution of the measurement of polaron motion.
The presented measurements are a compromise that optimizes the resolution and measurement area.Supplementary Fig. 10 shows the data of relaxation as a function of temperature for all measurements individually.Note that the scatter of the data is very small at low temperature, but increases with temperature.As discussed in the main text, the temperature-activated domain relaxation processes have a range of different   , reflecting the wide spread of relaxation rates, as defined by the changing energy landscape each time the configuration changes 2,3 .The average rate is presented in Fig. 1 of the main paper.The boundaries of the blue shaded area represent fits to the phenomenological model 13 given in the main text, where ) with values of   = 40  (top curve) and   = 200  (bottom curve) with  0 = 0.0023  −1 and 0.001  −1 respectively.The blue curve is the fit to the average of all the data shown in the main text, giving   = 36 , with  0 = 0.0011 −1 .Note that the lowtemperature value for the quantum rate   = 1.9 × 10 −4  −1 is the same in all three fits.The patterned area represents the base quantum rate, which is temperature independent.

Tip current and voltage dependence of domain reconfiguration: measurement details
The thermal calculation in section 1.2 has shown that during scanning the effect of Joule heating on the sample is negligible.However, the small tunneling currents may transfer sufficient charge to modify the charge configuration state according to the established mechanisms 1,14,15 .
To determine this empirically and determine the role of the tip in the relaxation process, we performed systematic measurements with different scanning parameters (I and V) both at 5 K and 50 K.
Exploring first the effect of STM bias voltage and tip-tunnelling current on R, we find that increasing either the tip current, or the tip-to-sample voltage at 5 K results in a slight increase on the relaxation rate (Supplementary Fig. 12a).Supplementary Fig. 11 shows an example of relaxation series at 5 K for 40-fold increased tunneling current 2 nA and -0.8 V, compared with 50 pA and -0.8 V for the sequences presented in Fig. 1 of the main text.The 40-fold increase in tip current (from 50 pA to 2 nA) increases R rather modestly, from ~2.5×10 -4 s -1 to ~4.2×10 -4 s - 1 (a factor of ~1.7).Similarly, quadrupling the tip voltage from -0.8 V to -3.2 V and simultaneously increasing the current three-fold from 50 pA to 150 pA results in a similarly modest change of R from ~2.5×10 -4 s -1 to ~5.5×10 -4 s -1 (a factor of 2.2).Supplementary Fig. 12 shows the comparison between the range of all rates measured for normal scanning parameters (50 pA and -0.8 V) and the average rates for two combinations of different parameters (2 nA and -0.8 V; 150 pA and -3.2 V).Thus, the application of huge increases in current or voltage apparently causes very small changes in R at 5 K.
At 50 K, the effect of changing the scanning parameters over the same range becomes undetectable.Supplementary Fig. 12b shows that relaxation rates at 50 K for different parameters (150 pA and -3.2 V; 2.0 nA and -0.8 V) are compatible with the rate at normal scanning parameters (50 pA and -0.8 V) presented in Fig. 1 of the main text.This confirms that the tip-induced relaxation is small in comparison with intrinsic quantum and thermally activated relaxation rates   and  0 respectively.
A possibility exists that external electrical noise influences R. The electrical noise floor within the UHV LT STM is many orders of magnitude below the value required to change the rate, as discussed above, so it is unlikely to cause enhanced relaxation.Cosmic rays are known to cause perturbations in quantum electrical devices, such as detectors and quantum processors (see main text), and they are also likely to trigger domain reconfiguration by locally depositing energy.However, this would be visible in the STM scans in the form of spikes.These are not very common on our apparatus, and cannot be correlated with relaxation, so at present cannot be identified as an external stimulus that drives relaxation.The presented data lend support to the arguments presented in the main text (and in the next section, 1.6) based on tunnelling arising from retracted-tip measurements (Fig. 1f).We emphasize that for the purpose of this study it is not important to eliminate completely external relaxation triggers, but only to ascertain that a finite intrinsic rate of domain reconfiguration is present.From our experiments presented above we can safely conclude that while the tip field or current and other external processes may enhance relaxation, a significant part of   arises from quantum processes that are intrinsic, and independent of external stimuli.
discussed in section 1.3) at the same scanning conditions (5 K, -0.8 V for tip-sample voltage bias, 50 pA for tunneling current).A histogram showing the polaron fraction moved in a number of repeated measurements is shown in Supplementary Fig. 15.
In order to estimate the tip-induced effects on the obtained results, the average fraction of moved polarons during dark relaxation was compared to the one measured during regular observation.The comparison shows that: This suggests that the fraction of polaron moved during the 8 + 32 minutes of dark relaxation can mostly be attributed to non-tip-induced processes.We conclude that domain reconfiguration takes place at a similar rate in the absence of tip scanning as with tip scanning.The tip thus has a relatively minor effect.
A fit to the 'dark' reconfiguration rate is shown in Fig. 1f of the main text, on the same timeline as for 'light' measurements where a scan takes place continuously.The rates are comparable.
2 Supplementary Notes 2 The low-temperature charge ordering in 1T-TaS2 is considered to be in the lattice Wigner crystal (WC) limit, which occurs when the electronic Coulomb energy significantly exceeds their kinetic energy 16 .The criterion for Wigner crystallization is defined in terms of the dimensionless Wigner-Seitz radius   = /, where  is the Coulomb energy and  is the kinetic energy.For a 2D electron system,   =  2 /(ℏ 2  1/2 ) , where  is the electron density,  is elementary charge, and  is the effective mass of electron.In 1T-TaS2,   is enhanced by the small carrier concentration (one electron per 13 unit cells -see Fig. 1a) and due to strong coupling to the lattice, polaronic effects renormalize the electronic bands, resulting in an increase of the effective mass.Altogether, this leads to a large expected   , and WC physics 16 .Doping, external charge injection 4,5,7,8 or optical perturbation 1,5 introduces additional charges in the WC, which leads to a complex phase diagram with the formation of domains, discommensurations 5 and charge fractionalization 14,16 (Fig. 1b).The resulting domain wall structure allows multiple possible configurations, arising from different relative domain positions, each of which is a well-defined macroscopic quantum state in the occupation number basis, separated from each other by a potential energy barrier.Such DWs arising from electron "overcrowding" have been observed experimentally by STM in numerous transition metal dichalcogenides 16 .1T-TaS2 is of particular interest because the domain reconfiguration dynamics can be tuned by temperature 5 such that at sufficiently low temperatures, the domain reconfigurations are slow enough to be recorded by scanning tunneling microscopy (STM).

Details of the theoretical model
In this chapter we introduce our theoretical model of 1T-TaS2 already developed and explored in our previous work 16 .We present its classical charge ordering aspects and their relevance to the simulation on D-Wave's quantum computer.The Hamiltonian presented in this section is mapped onto the D-Wave Hamiltonian   .
At the Fermi level, 1T-TaS2 exhibits a sparsely populated ( + ∑  q ( q †  q + 1/2) q , where  , is the hopping integral between lattice sites  and ,  the chemical potential,  , the annihilation operator of an electron at site  with spin ,  , =  , †  , .Electrons interact with each other via Coulomb interaction   (, ), as well as with phononic degrees of freedom, where  q is the annihilation operator of a phonon with wave vector q and  q is it's frequency.The electron-phonon interaction is taken to be Fröhlich like and it's matrix element is  , (q) and ℏ = 1.In the strong electron-phonon interaction limit,  can be solved exactly via the Lang Firsov canonical transformation ( ̃=    − ,  = ∑  , ( , (q) q − ℎ. . ) q,, ), where electrons get dressed by a surrounding lattice deformation and therefore become quasi particles dubbed as polarons.
We rewrite  ̃ with the classical 2D interaction Hamiltonian where (, ) =  0 exp(− , /  ) / , ,   is the screening radius,  , = |r  − r  |,   is the -th out of  lattice sites,   ∈ {0,1} is the occupation number of lattice site  and the sum runs over all lattice sites.The number of polarons in the system is varied by the chemical potential .
Polarons are considered as screened point charges that can only reside on Ta sites of the triangular atomic lattice.In this work we considered a triangular lattice with 2008 sites and open boundary conditions, where polarons repel each other via nearest neighbor repulsion. was set to a fixed value for which the ground state of   is a 1/3 polaronic lattice 16 .
describes the domain structure and salient features of the phase diagram very well 16,18  to the second term of   .The   term in   describes tunneling between different configurational states and with this additional term our simulated model can also be viewed as an extended two-dimensional TFIM on a triangular lattice.The triangular lattice with 2008 sites combined with the correlated structure of   requires a specific embedding which uses 2673 qubits subject to the TFIM.Details on the embedding can be found in a subsequent chapter.
In order to make the connection of our model as deployed on D-Wave and a Wigner crystal, we interpret the   term as the kinetic energy of polarons.Using the 2D Jordan Wigner transformation In the limit of strong interactions, domain walls dominate the system and the kinetic energy term simplifies into a hopping term for the domain walls ∑     †   + ℎ. .
The Wigner-Seitz radius   in our case is therefore related to the ratio between the second and the first term.The kinetic energy has an anomalous form and can be interpreted as the creation or annihilation of a spinless polaron at site , which depends on the current polaronic configuration.The comparison between the effect of this form of kinetic energy and the conventional hopping form is beyond the scope of this work and remains an interesting question for future studies.
In this paper we use an analogue to the 1/13 electronic lattice which forms in TAS.We use a 1/3 electronic lattice due to the smaller system sizes considered here.We have shown already in 16 that classically the system behaves very similar to the 1/13 case albeit at a higher polaron density.We assume here that this generalizes also to the quantum version of this model, which is considered in this work.
In order to capture non-equilibrium quantum domain wall reconfiguration dynamics, we require at least a few hundred atomic sites or qubits, which is impossible with state-of-the-art classical simulation methods and hardware 20 .Another promising candidate are quantum Monte Carlo simulations with added noise, which are not developed for directly probing non-equilibrium dynamics and it is still an open question whether they could be applied in our case.We also checked whether classical Monte Carlo simulations (next section) in the presence of a thermal bath can predict a crossover in ().As expected, thermally activated processes are responsible for the -dependent reconfiguration rate.However, they cannot explain the -independent behavior we observe in 1T-TaS2 at low  for interaction potentials known from previous literature 14,21,22 .

Classical Monte Carlo simulations
We performed classical Monte Carlo (MC) simulations in the form of simulated annealing, analogous to the quantum annealing performed on the D-Wave Advantage.Our classical simulations consist of Markov chain dynamics with the Hamiltonian   = , on a triangular lattice using the standard Metropolis algorithm.The proposed successive state in the chain consists of switching the positions of a particle and a hole on two randomly selected nearest neighboring sites.After this step, we also select a random site and destroy a particle if it's occupied or create one if it's empty.The simulations were performed at a fixed value of  =  0 on a triangular lattice with 2008 sites and lattice spacing , depicted in the section "Applying the model on the D-Wave's machine".The interaction potential (, ) is a positive constant  0 when  and  are nearest neighbors and, in a subsequent case, also nextnearest neighbors.Simulated annealing was performed by decreasing the effective temperature   =   / 0 from 0.5 to 0.001 in a 100 steps with 1000 MC sweeps at each temperature point.Our system exhibits an ordering phase transition from a high-temperature disordered gaseous phase to a crystalline ordered phase at low temperature.Supplementary Fig. 16 shows snapshots of the simulated annealing process at 4 different values of   .The left-most image is far above the phase transition, where the polaronic configuration is essentially random.The next two images toward the right are close above and below the phase transition temperature.The right-most image shows the final state, which is a striped ordered state.This is not the ground state of the model, which is a uniform blue lattice.The reason why we don't observe the ground state is due to the fact that the domain walls cost zero energy in terms of interaction.However, the energy could be lowered still with the addition of more polarons in the domain walls and thereby decreasing the chemical potential term.When the simulation crosses the phase transition, bubbles of the 3 possible ground state lattices begin to form and inevitably produce domain walls.After their formation they are very hard to get rid of from the point of view of simulated annealing, due to their low energy cost.For example, the final state in Supplementary Fig. 16 has an energy per particle (if  0 = 1) −0.7316 compared to the −0.9895 in the case of the same Hamiltonian being applied to the right-most configuration in Supplementary Fig. 18.
Despite our faliure of finding the ground state using MC simulations, it is still important to investigate this specific instance of the Hamiltonian and observe if the rate of polaronic movements (  ) can explain the experimentally observed rate.Supplementary Fig. 17 shows the average (  ) during the MC simulation.There is an apparent phase transition at   ≈ 0.19 when the slope of the rate changes abruptly.Also, the rate saturates above the phase transition due to the system effectively sampling from random polaronic configurations.It saturates at low   due to the system's finding a local minimum from where it cannot escape.There are are a few similarities and key differences between the MC simulations, D-Wave simulations and the 1T-TaS2 experiment: () There is a saturation above the phase transition in both the MC and D-Wave simulations due to melting of the ordered phase and thereby sampling from effectively random configurations.There is no such saturation present in 1T-TaS2 because the crystalline order is never melted there.polaronic dynamics a priori.However, we can compare it with D-Wave simulations.We obtain them from plots like the one shown in Fig. 2e shown in main text and Supplementary Fig. 33.The dimensionless gaps are  =1 = 3.5 ± 0.3;  =2 = 5.3 ± 0.9;  =5 = 7.6 ± 0.9;  =10 = 13.8 ± 0.5;  =15 = 17.2 ± 1, where the subscript denotes the value of  in the simulation.The value increases with  as we would expect when the kinetic energy of polarons is reduced compared to the potential energy.The corresponding MC slope can be found in the region  < 1, which is inconsistent from the extrapolated value of  = 58 in the main text.Therefore, we cannot compare the results obtained from D-Wave simulations with results from MC simulations, which makes sense as the former is a quantum version of the latter.
In conclusion, we have found that MC simulations are highly inconsistent with experimental observations.The main reason is in the simple fact that a classical rate generally saturates to the value 0 under a certain temperature, and the quantum rate doesn't.We have shown that there is in principle a way for the classical rate to be non-zero for any temperature, but this is only possible in the specific case when the energy cost for movement is exactly 0. This is inconsistent with previous work done on 1T-TaS2.We have also shown that the quantum model deployed on D-Wave cannot be reproduced with MC simulations because the phase transition characteristics do not match.

WKB Approximation and Macroscopic Quantum Tunneling
For coherent quantum tunneling, which we do not expect to observe in either system, the WKB rate, for   >    is given by: where  is the width of the barrier, and  * is the effective mass of the object.The physical meaning of the attempt frequency  is related to the tunneling of the entire object (domain) of mass  * from one configuration to another in which all electrons change position by ~0.3 . it may be estimated by the energy at which the oscillatory dynamics of a domain patch crosses over from thermally activated to quantum, and the thermal energy becomes comparable to the energy of a quantum oscillator   ~1 2  *  2  2 .The reconfiguration we observe may involve 100 − 1000 polarons, depending on the size of the domains.Taking an individual polaron mass of   ≃ 3   , we may have  * ~ 3 × 10 3   .For two domains to fuse, the electrons collectively move a distance of approximately elementary crystal one unit cell,  ≃ 0.3 .At With this values of , ssuming a barrier width of 1 , and an experimentally determined   = 2 , gives an estimated   ≃ 3 × 10 −3  −1 .The rate estimate is very sensitive to the values in the exponent.However, the reasonable agreement between the experimentally measured   = 2.2 ± 0.3 × 10 −4  −1 and the estimated   above shows that macroscopic quantum tunnelling between domain configurations is plausible on this timescale for our system with reasonable input parameters.In this estimate we ignore dislocations in the domain structure that topologically hinder relaxation.Such topological defects are an emergent property, which arise from the microscopic interactions that MQT does not describe.We also ignore the effects of pinning of electrons due to crystal lattice imperfections that may reduce the rate of relaxation.

Quantum Kramers theory of metastable decay
A natural model for domain reconfiguration dynamics in 1T-TaS2 is a succession of single escapes from local minima in the energy landscape towards the global minimum (ground state), approximated locally by a cubic potential, assuming also some degree of coupling to an external Ohmic bath.Since the experimental Γ() is averaged over many different escape events at temperature , we may apply the canonical single particle escape from a cubic potential with Ohmic dissipation in order to estimate its average parameters.This relatively simple model predicts the classical high temperature Arrhenius law as well as the crossover to a temperature independent quantum decay at some crossover temperature  0 .The former arises from classical thermally activated hopping over the potential barrier   and the latter from tunneling events in the presence of dissipation, which exponentially suppresses the tunneling probability.From the experimental crossover temperature  0 = 20  and crossover temperature range Δ = 10 , we obtain according to Ref. 23 the dimensionless Ohmic damping parameter  = 0.3 and the ratio of the activation barrier energy over the curvature (energy scale) of the local metastable minimum   ℏ 0 ≃ 0.4 ∼ 1.4, which is inconsistent with the underlying assumption of this model   /ℏ 0 ≫ 1.Therefore, quantum Kramers theory of metastable decay cannot explain the relaxation dynamics observed in 1T-TaS2.

Details on the incoherent macroscopic tunneling of our model
Here we make the argument for the mechanism of the decay process in 1T-TaS2 being incoherent macroscopic tunneling.We observe in experiment that the system starts off in some non-equilibrium domain state.Afterwards, we observe a reconfiguration event on a timescale of ∼ 10 3 , which is separated by a factor of 10 15 from the timescale obtained from the Arrhenius fit to the experimental rate ℏ/  ∼ 1 .It follows that the system is stuck in the initial state for a long time before decaying to the next configuration.Therefore, only the local energy levels in the energy spectrum are relevant to the decay, the next level lower in energy from the initial level being the most important.This is why a local two-level system approximation is appropriate for each observed reconfiguration event.The whole relaxation process can thereby be thought of as a succession of local decays in a two-level system, which is characterized by the tunneling rate Γ  ∼ Δ 2 / (discussed below).This is shown schematically in Fig. 2f.First, we assume a two level system with a bare gap Δ.If we couple such a system to a two-level oscillator bath, the 2 levels will split into 4 and the gap is reduced by ∼ .If the number of oscillators tends to ∞, the 2 levels become bands with width .In our consideration of both the quantum annealer, as well as 1T-TaS2, we assume a large enough  for the two energy bands to overlap.Therefore, the transitions of the oscillators of the bath are the cause of the incoherent macroscopic tunneling processes we are observing.Now we investigate the consequences of the conjecture that incoherent macroscopic tunneling is the mechanism responsible for domain reconfiguration dynamics in 1T-TaS2, driven by the Hamiltonian  =   +   +   , where   describes the polaronic system,   represents a generic harmonic oscillator bath at temperature , described by the spectral density () and   the system-bath interaction.If a two-level system is chosen as an approximation for a single reconfiguration event from a higher to a lower energy level in the full polaronic spectrum, where   = −(  + Δ  )/2, where   are Pauli matrices,  is the bias towards the system being in one of the two levels represented by the 2 values of   and Δ the matrix element for pure tunneling between the two levels, analytical results are known from literature [24][25][26][27][28] .According to Amin et al 25 we assume   =  ⨂   , where  is described by the bath spectral density , where ⟨… ⟩ denotes averaging over the bath degrees of freedom, and write their main result Γ q ∼ Δ 2 /, where  2 = ∫  2 () is the bandwidth of the bath noise.Γ  is not the exact result of Amin et al; it only represents the scale of the quantum decay rate they obtained for level transitions.If the system is initialized in some nonequilibrium state, the bath will of course drive the system towards its thermal equilibrium at temperature , but the timescale of these thermally induced transitions is still determined by Γ  .Transitions occur between different system-bath product states |〉|〉 and are characterized by an incoherent tunneling rate Γ  ∼ Δ 2 /.For example, in a single transition |〉 might stay the same, while |〉 undergoes some change, which leads to dephasing or decoherence of the whole system-bath state.This means that it is entirely possible for |〉 states to retain coherence and for us to not observe any transitions, even though meanwhile, the bath is dephasing the whole state.On the other hand, |〉 might undergo a transition simply due to a transition in |〉.This occurs when the two bands corresponding to the two spin states are overlapping, which means that for some energy region, |〉|〉 states are separated by the energy gaps between |〉 states and not Δ.The rate for such transitions is Γ  .According to the aforementioned equivalence assumption that 1T-TaS2 (M) and the quantum annealer (P) are both governed by the same , and the fact that we are able to reproduce the saturation of () in M with P, meanwhile observing similar reconfiguration processes in P as in M, we argue that the perturbation theory developed for P is also applicable to M. Now we apply the developed perturbation theory to our experiment and simulations, which takes the system's energy landscape into account through considering the entire energy spectrum system and focusing on the transitions between different energy levels (Fig. 2f).We assume that when the STM tip measures a polaronic configuration, it measures the occupation number ⟨  |  |  ⟩ for every atomic site  and is insensitive to the bath degrees of freedom, which represent the noise of the system.The unperturbed Hamiltonian (when the STM tip is present) is therefore the potential energy between polarons, the coupling to the bath, and the bath (()  +   +   ), because they only contain   and bath operators.The unperturbed energy spectrum contains all the possible polaronic configurations ranked by potential energy, which are broadened into bands by the bath.Of course, a single STM measurement forces the system into one specific polaronic configuration, where the bath degrees of freedom remain unspecified.This means that the whole polaron-bath system can lie anywhere on one of the polaronic energy spectrum bands, where the bath degrees of freedom change across the band, but the polaronic configuration is the same.The perturbation is the kinetic energy of polarons (()  ), which is always present but only contributes significantly when the STM is removed or turned off.Since the time-scale   for relaxation is so large compared to   , the probability for transitioning to another configurational state is ≪ 1, thus justifying the perturbation theory approach.
The gap Δ = () is characterized by the overall energy scale of the system , where () is now the dimensionless gap, which typically depends exponentially on the Hamiltonian parameters 29,30 .The energy scale of P,   = 1.6 ± 0.8 , governing the low temperature relaxation dynamics of system P is set by the magnetic fields operating on the SQUIDs serving as qubits on the QPU 31 , while   = 3 ± 0.4  was measured in the STM experiment.We also know from previous literature 32,33 that the spectral density of noise in both M and P is well described as 1/ noise at low frequencies, which we assume is the dominant effect of noise in both cases.In order to make the connection from the fundamental time-scale to quantum incoherent relaxation, we rewrite Γ  = .If we combine the same 1/ shape of the noise spectrum in both P and M with a similar quantitative value of energy scale to noise bandwidth ratio, we can conclude that the influence of noise is the same in P and M.
3 Supplementary Notes 3 In order to simulate quantum melting of domains we need to first be able to map our model on to an actual D-Wave machine.In this paper we use the newly released Advantage6.1 quantum computer with 5436 physical qubits.In this chapter we present all the necessary steps required for a successful deployment of our model.), (11)

1. 2 . 1
Thermal behavior of the sample during the creation of the domain state by a single electrical pulse Supplementary Figure 2 | Temperature (a) and current density (b) maps for the strong pulse used to create the domain state.The starting state parameters correspond to the high-resistivity C state.The domain state is examined at ~50 nm from the center.The current density is in units of Amperes/m 3 .

. 10 ,
of the energy scale, as well as the measured time-scale   ≈ 1000  ≈ 10 15  0  , where  0 .If Γ = 1/, then Γ  = 10 −15   /ℏ ∼ 10 −15   .The 15 orders of magnitude difference is therefore explained by the quantity () 2     = 10 −15 .In order to estimate the ratio     in system M, we turn to the definition  2 = ∫  2()27 and turn to measurements of low-frequency voltage fluctuations32 , which suggest that we can approximate the spectral density of the external bath of M with 1/ behavior as   (Integrating in the range of 3 − 100  returns a value of   ≈ 0.3 .Therefore, which leads to () ∼ 10 −8 , which in turn leads to a gap Δ  ∼ 10  between different polaronic configurational states.If we take the value of   from literature33 and calculate
, but not quantum dynamics.To introduce quantum tunneling between different configurations of domains we first simplify   to   ′ = ∑  ,      , , where  , = 1 2 (, ) −  , .This makes  ̃ applicable to the quadratic unconstrained binary optimization formalism implemented on the D-